%read frequency in gates file from flowjo tables
%Array: [row,col,timepoint,stat]
%stats = [Q1,Q2,Q3,Q4,Geo-FSC,Geo-SSC]
%flow_path = '../All_Lim_FACS_tvs/120512_freq_gates_96_from_111912.tab';
flow_path = '../All_Lim_FACS_tvs/121012_96all.tab';

[data,time,titles] = read_flowjo_table_time(flow_path);

%see input experiments
titles'

%Dox (ug/mL)
dox_tvs = [12.500	6.250	3.125	1.563	0.781	0.391	0.195	0.098	0.049	0.024	0.012	0.000];
dox_tvs = fliplr(dox_tvs);
dox_hyun = [0, 0.01, 0.02, 0.04, 0.08, 0.2, 0.4, 0.8, 2, 4, 8, 16];
%visualize dose axis differences
%{
figure, hold on
scatter(1:numel(dox_hyun),log10(dox_hyun),'b')
scatter(1:numel(dox_tvs),log10(dox_tvs),'r')
%OD axis
%}
od = [0 1/8 1/4 1/2 1 2 4 8];

%Normalize percentages to exclude Q4, convert to fractions
for i = 1:size(data,1)
    for j = 1:size(data,2)
        for k = 1:size(data,3)
                data(i,j,k,1:3) = data(i,j,k,1:3)/sum(data(i,j,k,1:3));
                data(i,j,k,4) = data(i,j,k,4)/100;
        end
    end
end

%Calculate max cluster in any experiment (for cmap)
clusters = data(:,:,:,2);
cluster_max = max(clusters(:));
%% Plotting
defaultPlot()

%Delete unused NaN vector space
i=1;
while ~all(all(all((isnan(squeeze(data(:,:,i,:)))))))
    i=i+1;
end
data(:,:,i:end,:)=[];

%Flip row/col of thomas's experiments to match Hyun's
%Flip row/col of a single Hyun experiment -- visualized error
    % no clusters in this experiment, removed
for i=1:numel(titles)
    if any(strfind(titles{i},'yTS'))
        for j=1:size(data,4)
            data(:,:,i,j) = flipdim(flipdim(data(:,:,i,j),1),2);
        end
    end
    if any(strfind(titles{i},'BookIII 27 t19'))
        for j=1:size(data,4)
            data(:,:,i,j) = flipdim(data(:,:,i,j),1);
        end
    end
end

%Exptrapolate NaNs
for m = 1:size(data,3)
    for n=1:size(data,4)
        datacopy = squeeze(data(:,:,m,n));
        [x,y] = meshgrid(1:12,1:8);
        x(isnan(datacopy))=[];
        y(isnan(datacopy))=[];
        datacopy(isnan(datacopy))=[];
        [xi,yi] = meshgrid(1:12,1:8);
        F = TriScatteredInterp(x(:),y(:),datacopy(:),'nearest');
        data(:,:,m,n) = F(xi,yi);
    end
end

%Define White->Dark Blue colormap
blue_clim = [[1;1;1] [0;0;1] [0;0;0]];
blue_c = [0 128 256];
blue_ci = [1:256];
for i=1:length(blue_c)
    blue_cmap(i,:) = interp1(blue_c,blue_clim(i,:),blue_ci); 
end
blue_cmap=blue_cmap';

%Check all contour maps for orientation
%for i=1:numel(titles)
for i = 5:9
    figure
    contourf(dox_tvs,od,squeeze(data(:,:,i,2)),'LineColor','w')
    colormap(blue_cmap)
    %xlabel('Dox...')
    set(gca(),'XScale','log')
    xlim([1e-3 12.5])
    set(gca(),'XTick',[1e-3 1e-1 1 1e1])
    set(gca(),'XTickLabel',{'(No dox)',1e-1,1,1e1})
    %ylabel('OD Ratio')
    set(gca(),'YScale','log')
    ylim([1e-2 8])
    set(gca(),'YTick',[1e-2 .125 1 8])
    set(gca(),'YTickLabel',{'(Red only)','1 : 8','1 : 1','8 : 1'})
    title(titles{i}(2:end));
    caxis([0 cluster_max])
    colorbar()
    set(gca,'FontSize',24)
end
%}

%separate independent experiments
%Trajectory = [initial,vector index,timept,experiment]
traj = nan(size(data,1), 3, max(ends-starts+1), length(starts));
starts = [1 find(diff(time)<0)+1];
ends = [find(diff(time)<0) length(time)];
for ex = 1:length(starts)
    for t = 1:ends(ex)-starts(ex)+1
        traj(:,:,t,ex) = data(:,end,starts(ex)+t-1,1:3);
    end
end

%Need to replace clusters->3, green->2 (currently r-c-g)
temp = traj(:,2,:,:);
traj(:,2,:,:) = traj(:,3,:,:);
traj(:,3,:,:) = temp;

%Check trajectories in 3D
for i=1:size(traj,4)
    figure, hold on
    %plot3(squeeze(traj(:,1,:,i))',squeeze(traj(:,2,:,i))',squeeze(traj(:,3,:,i))');
end
%set as Traj = [vector,timept,initial,experiment]
traj = permute(traj,[2 3 1 4]);
%Transformation matrix onto 2-simplex
T = [0 1 .5;0 0 1];

for ic = 1:size(traj,3)
    for ex = 1:size(traj,4)
        trajT(:,:,ic,ex) = T*traj(:,:,ic,ex);
    end
end
cmap = [fliplr(linspace(0,1,8));linspace(0,1,8);[0:1/3:1 1:-1/3:0]]';
figure
for ex = [1 2 4]
    hold on
    for ic = 1:size(traj,3) 
        
        %plot(trajT(1,:,ic,ex),trajT(2,:,ic,ex))
        x = trajT(1,:,ic,ex);
        y = trajT(2,:,ic,ex);
        quiver(x,y,[diff(x) 0],[diff(y) 0],0,'LineWidth',2,'Color',cmap(ic,:))
        axis([0 1 0 1])
    end
end
%draw simplex
line([0 1],[0 0],'Color','k','LineWidth',2)
line([1 .5],[0 1],'Color','k','LineWidth',2)
line([.5 0],[1,0],'Color','k','LineWidth',2)